!> \file ysuvdif.F90
!! This file contains the CCPP-compliant YSU scheme which computes
!! subgrid vertical turbulence mixing using traditional K-profile method
!! Please refer to (Hong, Noh and Dudhia, 2006, MWR).
!!
!! Subroutine 'ysuvdif_run' computes subgrid vertical turbulence mixing
!! using YSU K-profile method
!!
!----------------------------------------------------------------------

      module ysuvdif
      contains

      subroutine ysuvdif_init (do_ysu,errmsg,errflg)

        logical,          intent(in) :: do_ysu
        character(len=*), intent(out) :: errmsg
        integer,          intent(out) :: errflg

        ! Initialize CCPP error handling variables
        errmsg = ''
        errflg = 0

        ! Consistency checks
        if (.not. do_ysu) then
          write(errmsg,fmt='(*(a))') 'Logic error: do_ysu = .false.'      
          errflg = 1
          return
        end if
      end subroutine ysuvdif_init

!> \defgroup YSU FV3GFS ysuvdif_run Main
!! \brief This subroutine contains all of the logic for the
!! YSU scheme.
!!
!> \section arg_table_ysuvdif_run Argument Table
!! \htmlinclude ysuvdif_run.html
!!
!-------------------------------------------------------------------------------
   subroutine ysuvdif_run(im,km,ux,vx,tx,qx,p2d,p2di,pi2d,karman,              &
                    utnp,vtnp,ttnp,qtnp,                                       &
                    swh,hlw,xmu,ntrac,ndiff,ntcw,ntiw,                         &
                    phii,phil,psfcpa,                                          &
                    zorl,stress,hpbl,psim,psih,                                &
                    landmask,heat,evap,wspd,br,                                &
                    g,rd,cp,rv,ep1,ep2,xlv,                                    &
                    dusfc,dvsfc,dtsfc,dqsfc,                                   &
                    dt,kpbl1d,u10,v10,lssav,ldiag3d,qdiag3d,                   &
                    flag_for_pbl_generic_tend,ntoz,ntqv,dtend,dtidx,           &
                    index_of_temperature,index_of_x_wind,index_of_y_wind,      &
                    index_of_process_pbl,errmsg,errflg   )

   use machine , only : kind_phys
!
!-------------------------------------------------------------------------------
   implicit none
!-------------------------------------------------------------------------------
   real(kind=kind_phys),parameter    ::  xkzminm = 0.1,xkzminh = 0.01
   real(kind=kind_phys),parameter    ::  xkzmin = 0.01,xkzmax = 1000.,rimin = -100.
   real(kind=kind_phys),parameter    ::  rlam = 30.,prmin = 0.25,prmax = 4.
   real(kind=kind_phys),parameter    ::  brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
   real(kind=kind_phys),parameter    ::  afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0
   real(kind=kind_phys),parameter    ::  phifac = 8.,sfcfrac = 0.1
   real(kind=kind_phys),parameter    ::  d1 = 0.02, d2 = 0.05, d3 = 0.001
   real(kind=kind_phys),parameter    ::  h1 = 0.33333335, h2 = 0.6666667
   real(kind=kind_phys),parameter    ::  zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
   real(kind=kind_phys),parameter    ::  tmin=1.e-2
   real(kind=kind_phys),parameter    ::  gamcrt = 3.,gamcrq = 2.e-3
   real(kind=kind_phys),parameter    ::  xka = 2.4e-5
   real(kind=kind_phys),parameter    ::  rcl = 1.0
   real(kind=kind_phys),intent(in)   ::  karman
   integer,parameter ::  imvdif = 1
   integer,parameter ::  ysu_topdown_pblmix = 1
!
!-------------------------------------------------------------------------------------
!  input variables
   integer,  intent(in   )   ::     im,km,ntrac,ndiff,ntcw,ntiw,ntoz
   real(kind=kind_phys),     intent(in   )   ::     g,cp,rd,rv,ep1,ep2,xlv,dt

   real(kind=kind_phys),     dimension( :,: ),                                    &
             intent(in)      ::                 pi2d,p2d,phil,ux,vx,swh,hlw,tx

   real(kind=kind_phys),     dimension( :,:,: )                             , &
             intent(in   )   ::                                             qx

   real(kind=kind_phys),     dimension( :,: )                                , &
             intent(in   )   ::                                      p2di,phii

   real(kind=kind_phys),     dimension( : )                                      , &
             intent(in)  ::     stress,zorl,heat,evap,wspd,br,psim,psih,psfcpa,     &
                                                                   u10,v10,xmu
   integer,  dimension(:)                                                         ,&
             intent(in   )   ::                                      landmask
   logical,  intent(in   )   :: lssav, ldiag3d, qdiag3d,                            &
                                flag_for_pbl_generic_tend
!
!----------------------------------------------------------------------------------
! input/output variables
!
   real(kind=kind_phys),     dimension( :,: )                                   , &
             intent(inout)   ::                                utnp,vtnp,ttnp
   real(kind=kind_phys),     dimension( :,:,: )                              , &
             intent(inout)   ::                                          qtnp
   real(kind=kind_phys), optional, intent(inout) :: dtend(:,:,:)
   integer, intent(in) :: dtidx(:,:), ntqv, index_of_temperature,                  &
        index_of_x_wind, index_of_y_wind, index_of_process_pbl
!
!---------------------------------------------------------------------------------
! output variables
   integer,  dimension( : ), intent(out  )   ::                       kpbl1d
   real(kind=kind_phys),     dimension( : ),                                   &
             intent(out)   ::                                            hpbl
   real(kind=kind_phys),    dimension( : ),                                    &
              intent(out)  :: dusfc,dvsfc, dtsfc,dqsfc

   ! error messages
   character(len=*), intent(out)    ::                                 errmsg
   integer,          intent(out)    ::                                 errflg
!
!--------------------------------------------------------------------------------
!
! local vars
!
   real(kind=kind_phys),     dimension( im )            ::                hol
   real(kind=kind_phys),     dimension( im, km+1 ) ::                      zq
!
   real(kind=kind_phys),     dimension( im, km )   ::                          &
                                                               thx,thvx,thlix, &
                                                                          del, &
                                                                          dza, &
                                                                          dzq, &
                                                                        xkzom, &
                                                                        xkzoh, &
                                                                           za
!
   real(kind=kind_phys),    dimension( im )             ::                     &
                                                                         rhox, &
                                                                       govrth, &
                                                                  zl1,thermal, &
                                                                       wscale, &
                                                                  hgamt,hgamq, &
                                                                    brdn,brup, &
                                                                    phim,phih, &
                                                                        prpbl, &
                                                              wspd1,thermalli
!
   real(kind=kind_phys),    dimension( im, km )    ::                     xkzm,xkzh, &
                                                                        f1,f2, &
                                                                        r1,r2, &
                                                                        ad,au, &
                                                                           cu, &
                                                                           al, &
                                                                         xkzq, &
                                                                         zfac, &
                                                                        rhox2, &
                                                                       hgamt2
!
   real(kind=kind_phys),    dimension( im )    ::                                         &
                                                                         brcr, &
                                                                        sflux, &
                                                                         zol1, &
                                                                    brcr_sbro
!
   real(kind=kind_phys),    dimension( im )    ::                       xland
   real(kind=kind_phys),    dimension( im )    ::                         ust
   real(kind=kind_phys),    dimension( im )    ::                         hfx
   real(kind=kind_phys),    dimension( im )    ::                         qfx
   real(kind=kind_phys),    dimension( im )    ::                         znt
   real(kind=kind_phys),    dimension( im )    ::                         uox
   real(kind=kind_phys),    dimension( im )    ::                         vox
!
   real(kind=kind_phys),    dimension( im, km, ndiff)  ::                     r3,f3
   integer, dimension( im )             ::                  kpbl,kpblold
!
   logical, dimension( im )             ::                        pblflg, &
                                                                       sfcflg, &
                                                                       stable, &
                                                                     cloudflg

   logical                                   ::                     definebrup
!
   integer ::  n,i,k,l,ic,is,kk
   integer ::  klpbl, ktrace1, ktrace2, ktrace3
!
!
   real(kind=kind_phys)    ::  dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
   real(kind=kind_phys)    ::  ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
   real(kind=kind_phys)    ::  brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
   real(kind=kind_phys)    ::  utend,vtend,ttend,qtend
   real(kind=kind_phys)    ::  dtstep,govrthv
   real(kind=kind_phys)    ::  cont, conq, conw, conwrc, rovcp
!

   real(kind=kind_phys), dimension( im, km )     ::                wscalek,wscalek2
   real(kind=kind_phys), dimension( im )              ::                           wstar
   real(kind=kind_phys), dimension( im )              ::                           delta
   real(kind=kind_phys), dimension( im, km )     ::                     xkzml,xkzhl, &
                                                               zfacent,entfac
   real(kind=kind_phys), dimension( im )              ::                            ust3, &
                                                                       wstar3, &
                                                                     wstar3_2, &
                                                                  hgamu,hgamv, &
                                                                      wm2, we, &
                                                                       bfxpbl, &
                                                                hfxpbl,qfxpbl, &
                                                                ufxpbl,vfxpbl, &
                                                                        dthvx
   real(kind=kind_phys)    ::  prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx,                           &
               dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr,     &
               prfac,prfac2,phim8z,radsum,tmp1,templ,rvls,temps,ent_eff,    &
               rcldb,bruptmp,radflux
   integer                 ::  idtend
!
!-------------------------------------------------------------------------------
!
! Initialize CCPP error handling variables
   errmsg = ''
   errflg = 0

   klpbl = km
!
   rovcp=rd/cp
   cont=cp/g
   conq=xlv/g
   conw=1./g
   conwrc = conw*sqrt(rcl)
   conpr = bfac*karman*sfcfrac
!
!  change xland values
   do i=1,im
     if(landmask(i).eq.0) then !ocean
       xland(i) = 2
     else
       xland(i) = 1  !land
     end if
   end do
!
   do k = 1,km
     do i = 1,im
       thx(i,k) = tx(i,k)/pi2d(i,k)
       thlix(i,k) = (tx(i,k)-xlv*qx(i,k,ntcw)/cp-2.834E6*qx(i,k,ntiw)/cp)/pi2d(i,k)
     enddo
   enddo
!
   do k = 1,km
     do i = 1,im
       tvcon = (1.+ep1*qx(i,k,1))
       thvx(i,k) = thx(i,k)*tvcon
     enddo
   enddo
!
   do i = 1,im
     tvcon = (1.+ep1*qx(i,1,1))
     rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon)
     govrth(i) = g/thx(i,1)
     hfx(i) = heat(i)*rhox(i)*cp ! reset to the variable in WRF
     qfx(i) = evap(i)*rhox(i)    ! reset to the variable in WRF
     ust(i) = sqrt(stress(i))    ! reset to the variable in WRF
     znt(i) = 0.01*zorl(i)       ! reset to the variable in WRF
     uox(i) = 0.0
     vox(i) = 0.0
   enddo
!
!-----compute the height of full- and half-sigma levels above ground
!     level, and the layer thicknesses.
!
   do i = 1,im
     zq(i,1) = 0.
   enddo
!
   do k = 1,km
     do i = 1,im
       zq(i,k+1) = phii(i,k+1)*conw
       tvcon = (1.+ep1*qx(i,k,1))
       rhox2(i,k) = p2d(i,k)/(rd*tx(i,k)*tvcon)
     enddo
   enddo
!
   do k = 1,km
     do i = 1,im
       za(i,k) = phil(i,k)*conw
       dzq(i,k) = zq(i,k+1)-zq(i,k)
       del(i,k) = p2di(i,k)-p2di(i,k+1)
     enddo
   enddo
!
   do i = 1,im
     dza(i,1) = za(i,1)
   enddo
!
   do k = 2,km
     do i = 1,im
       dza(i,k) = za(i,k)-za(i,k-1)
     enddo
   enddo

!    write(0,*)"===CALLING ysu; input:"
!         print*,"t:",tx(1,1),tx(1,2),tx(1,km)
!         print*,"u:",ux(1,1),ux(1,2),ux(1,km)
!         print*,"v:",vx(1,1),vx(1,2),vx(1,km)
!         print*,"q:",qx(1,1,1),qx(1,2,1),qx(1,km,1)
!         print*,"exner:",pi2d(1,1),pi2d(1,2),pi2d(1,km)
!         print*,"phii:",zq(1,1),zq(1,2),zq(1,km+1)
!         print*,"phil:",za(1,1),za(1,2),za(1,km)
!         print*,"p2d:",p2d(1,1),p2d(1,2),p2d(1,km)
!         print*,"p2di:",p2di(1,1),p2di(1,2),p2di(1,km+1)
!         print *,"del:",del(1,1),del(1,2),del(1,km)
!         print*,"znt,ust,wspd:",znt(1),ust(1),wspd(1)
!         print*,"hfx,qfx,xland:",hfx(1),qfx(1),xland(1)
!         print*,"rd,rv,g:",rd,rv,g
!         print*,"ep1,ep2,xlv:",ep1,ep2,xlv
!         print*,"br,psim,psih:",br(1),psim(1),psih(1)
!         print*,"u10,v10:",u10(1),v10(1)
!         print*,"psfcpa,cp:",psfcpa(1),cp
!         print*,"ntrac,ndiff,ntcw,ntiw:",ntrac,ndiff,ntcw,ntiw
!
!
!-----initialize vertical tendencies and
!
!   utnp(:,:) = 0.
!   vtnp(:,:) = 0.
!   ttnp(:,:) = 0.
!   qtnp(:,:,:) = 0.
!
   do i = 1,im
     wspd1(i) = sqrt( (ux(i,1)-uox(i))*(ux(i,1)-uox(i)) + (vx(i,1)-vox(i))*(vx(i,1)-vox(i)) )+1.e-9
   enddo
!
!---- compute vertical diffusion
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     compute preliminary variables
!
   dtstep = dt
   dt2 = 2.*dtstep
   rdt = 1./dt2
!
   do i = 1,im
     bfxpbl(i) = 0.0
     hfxpbl(i) = 0.0
     qfxpbl(i) = 0.0
     ufxpbl(i) = 0.0
     vfxpbl(i) = 0.0
     hgamu(i)  = 0.0
     hgamv(i)  = 0.0
     delta(i)  = 0.0
     wstar3_2(i) =  0.0
   enddo
!
   do k = 1,klpbl
     do i = 1,im
       wscalek(i,k) = 0.0
       wscalek2(i,k) = 0.0
     enddo
   enddo
!
   do k = 1,klpbl
     do i = 1,im
       zfac(i,k) = 0.0
     enddo
   enddo
   do k = 1,klpbl-1
     do i = 1,im
       xkzom(i,k) = xkzminm
       xkzoh(i,k) = xkzminh
     enddo
   enddo
!
   do i = 1,im
     dusfc(i) = 0.
     dvsfc(i) = 0.
     dtsfc(i) = 0.
     dqsfc(i) = 0.
   enddo
!
   do i = 1,im
     hgamt(i)  = 0.
     hgamq(i)  = 0.
     wscale(i) = 0.
     kpbl(i)   = 1
     hpbl(i)   = zq(i,1)
     zl1(i)    = za(i,1)
     thermal(i)= thvx(i,1)
     thermalli(i) = thlix(i,1)
     pblflg(i) = .true.
     sfcflg(i) = .true.
     sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1)
     if(br(i).gt.0.0) sfcflg(i) = .false.
   enddo
!
!     compute the first guess of pbl height
!
   do i = 1,im
     stable(i) = .false.
     brup(i) = br(i)
     brcr(i) = brcr_ub
   enddo
!
   do k = 2,klpbl
     do i = 1,im
       if(.not.stable(i))then
         brdn(i) = brup(i)
         spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
         brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
         kpbl(i) = k
         stable(i) = brup(i).gt.brcr(i)
       endif
     enddo
   enddo
!
   do i = 1,im
     k = kpbl(i)
     if(brdn(i).ge.brcr(i))then
       brint = 0.
     elseif(brup(i).le.brcr(i))then
       brint = 1.
     else
       brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
     endif
     hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
     if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
     if(kpbl(i).le.1) pblflg(i) = .false.
   enddo
!
   do i = 1,im
     fm = psim(i)
     fh = psih(i)
     zol1(i) = max(br(i)*fm*fm/fh,rimin)
     if(sfcflg(i))then
       zol1(i) = min(zol1(i),-zfmin)
     else
       zol1(i) = max(zol1(i),zfmin)
     endif
     hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac
     if(sfcflg(i))then
       phim(i) = (1.-aphi16*hol1)**(-1./4.)
       phih(i) = (1.-aphi16*hol1)**(-1./2.)
       bfx0 = max(sflux(i),0.)
       hfx0 = max(hfx(i)/rhox(i)/cp,0.)
       qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
       wstar3(i) = (govrth(i)*bfx0*hpbl(i))
       wstar(i) = (wstar3(i))**h1
     else
       phim(i) = (1.+aphi5*hol1)
       phih(i) = phim(i)
       wstar(i)  = 0.
       wstar3(i) = 0.
     endif
     ust3(i)   = ust(i)**3.
     wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
     wscale(i) = min(wscale(i),ust(i)*aphi16)
     wscale(i) = max(wscale(i),ust(i)/aphi5)
   enddo
!
!     compute the surface variables for pbl height estimation
!     under unstable conditions
!
   do i = 1,im
     if(sfcflg(i).and.sflux(i).gt.0.0)then
       gamfac   = bfac/rhox(i)/wscale(i)
       hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
       hgamq(i) = min(gamfac*qfx(i),gamcrq)
       vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
       thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
       thermalli(i)= thermalli(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
       hgamt(i) = max(hgamt(i),0.0)
       hgamq(i) = max(hgamq(i),0.0)
       brint    = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
       hgamu(i) = brint*ux(i,1)
       hgamv(i) = brint*vx(i,1)
     else
       pblflg(i) = .false.
     endif
   enddo
!
!     enhance the pbl height by considering the thermal
!
   do i = 1,im
     if(pblflg(i))then
       kpbl(i) = 1
       hpbl(i) = zq(i,1)
     endif
   enddo
!
   do i = 1,im
     if(pblflg(i))then
       stable(i) = .false.
       brup(i) = br(i)
       brcr(i) = brcr_ub
     endif
   enddo
!
   do k = 2,klpbl
     do i = 1,im
       if(.not.stable(i).and.pblflg(i))then
         brdn(i) = brup(i)
         spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
         brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
         kpbl(i) = k
         stable(i) = brup(i).gt.brcr(i)
       endif
     enddo
   enddo
!
!     enhance pbl by theta-li
!
   if (ysu_topdown_pblmix.eq.1)then
     do i = 1,im
        kpblold(i) = kpbl(i)
        definebrup=.false.
        do k = kpblold(i), km-1
           spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
           bruptmp = (thlix(i,k)-thermalli(i))*(g*za(i,k)/thlix(i,1))/spdk2
           stable(i) = bruptmp.ge.brcr(i)
           if (definebrup) then
           kpbl(i) = k
           brup(i) = bruptmp
           definebrup=.false.
           endif
           if (.not.stable(i)) then !overwrite brup brdn values
           brdn(i)=bruptmp
           definebrup=.true.
           pblflg(i)=.true.
           endif
        enddo
     enddo
   endif

   do i = 1,im
     if(pblflg(i)) then
       k = kpbl(i)
       if(brdn(i).ge.brcr(i))then
         brint = 0.
       elseif(brup(i).le.brcr(i))then
         brint = 1.
       else
         brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
       endif
       hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
       if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
       if(kpbl(i).le.1) pblflg(i) = .false.
     endif
   enddo
!
!     stable boundary layer
!
   do i = 1,im
     if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
       brup(i) = br(i)
       stable(i) = .false.
     else
       stable(i) = .true.
     endif
   enddo
!
   do i = 1,im
     if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then
       wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
       wspd10 = sqrt(wspd10)
       ross = wspd10 / (cori*znt(i))
       brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
     endif
   enddo
!
   do i = 1,im
     if(.not.stable(i))then
       if((xland(i)-1.5).ge.0)then
         brcr(i) = brcr_sbro(i)
       else
         brcr(i) = brcr_sb
       endif
     endif
   enddo
!
   do k = 2,klpbl
     do i = 1,im
       if(.not.stable(i))then
         brdn(i) = brup(i)
         spdk2   = max(ux(i,k)**2+vx(i,k)**2,1.)
         brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
         kpbl(i) = k
         stable(i) = brup(i).gt.brcr(i)
       endif
     enddo
   enddo
!
   do i = 1,im
     if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
       k = kpbl(i)
       if(brdn(i).ge.brcr(i))then
         brint = 0.
       elseif(brup(i).le.brcr(i))then
         brint = 1.
       else
         brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
       endif
       hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
       if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
       if(kpbl(i).le.1) pblflg(i) = .false.
     endif
   enddo
!
!     estimate the entrainment parameters
!
   do i = 1,im
     cloudflg(i)=.false. 
     if(pblflg(i)) then
       k = kpbl(i) - 1
       wm3       = wstar3(i) + 5. * ust3(i)
       wm2(i)    = wm3**h2
       bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
       dthvx(i)  = max(thvx(i,k+1)-thvx(i,k),tmin)
       we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
       if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3.and.ysu_topdown_pblmix.eq.1)then
           if ( kpbl(i) .ge. 2) then
                cloudflg(i)=.true. 
                templ=thlix(i,k)*(p2di(i,k+1)/100000)**rovcp
                !rvls is ws at full level
                rvls=100.*6.112*EXP(17.67*(templ-273.16)/(templ-29.65))*(ep2/p2di(i,k+1))
                temps=templ + ((qx(i,k,1)+qx(i,k,ntcw))-rvls)/(cp/xlv  + &
                ep2*xlv*rvls/(rd*templ**2))
                rvls=100.*6.112*EXP(17.67*(temps-273.15)/(temps-29.65))*(ep2/p2di(i,k+1))
                rcldb=max((qx(i,k,1)+qx(i,k,ntcw))-rvls,0.)
                !entrainment efficiency
                dthvx(i)  = (thlix(i,k+2)+thx(i,k+2)*ep1*(qx(i,k+2,1)+qx(i,k+2,ntcw))) &
                          - (thlix(i,k) + thx(i,k)  *ep1*(qx(i,k,1)  +qx(i,k,ntcw)))
                dthvx(i)  = max(dthvx(i),0.1)
                tmp1      = xlv/cp * rcldb/(pi2d(i,k)*dthvx(i))
                ent_eff   = 0.2 * 8. * tmp1 +0.2

                radsum=0.
                do kk = 1,kpbl(i)-1
                   radflux=swh(i,kk)*xmu(i)+hlw(i,kk) !radiative heating rate temp/s
                   radflux=radflux*cp/g*(p2di(i,kk)-p2di(i,kk+1)) ! converts temp/s to W/m^2
                   if (radflux < 0.0 ) radsum=abs(radflux)+radsum
                enddo
                radsum=max(radsum,0.0)

                !recompute entrainment from sfc thermals
                bfx0 = max(max(sflux(i),0.0)-radsum/rhox2(i,k)/cp,0.)
                bfx0 = max(sflux(i),0.0)
                wm3 = (govrth(i)*bfx0*hpbl(i))+5. * ust3(i)
                wm2(i)    = wm3**h2
                bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
                dthvx(i)  = max(thvx(i,k+1)-thvx(i,k),tmin)
                we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))

                !entrainment from PBL top thermals
                bfx0 = max(radsum/rhox2(i,k)/cp-max(sflux(i),0.0),0.)
                bfx0 = max(radsum/rhox2(i,k)/cp,0.)
                wm3       = (g/thvx(i,k)*bfx0*hpbl(i)) ! this is wstar3(i)
                wm2(i)    = wm2(i)+wm3**h2
                bfxpbl(i) = - ent_eff * bfx0
                dthvx(i)  = max(thvx(i,k+1)-thvx(i,k),0.1)
                we(i) = we(i) + max(bfxpbl(i)/dthvx(i),-sqrt(wm3**h2))

                !wstar3_2
                bfx0 = max(radsum/rhox2(i,k)/cp,0.)
                wstar3_2(i) =  (g/thvx(i,k)*bfx0*hpbl(i))
                !recompute hgamt 
                wscale(i) = (ust3(i)+phifac*karman*(wstar3(i)+wstar3_2(i))*0.5)**h1
                wscale(i) = min(wscale(i),ust(i)*aphi16)
                wscale(i) = max(wscale(i),ust(i)/aphi5)
                gamfac   = bfac/rhox(i)/wscale(i)
                hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
                hgamq(i) = min(gamfac*qfx(i),gamcrq)
                gamfac   = bfac/rhox2(i,k)/wscale(i)
                hgamt2(i,k) = min(gamfac*radsum/cp,gamcrt)
                hgamt(i) = max(hgamt(i),0.0) + max(hgamt2(i,k),0.0)
                brint    = -15.9*ust(i)*ust(i)/wspd(i)*(wstar3(i)+wstar3_2(i))/(wscale(i)**4.)
                hgamu(i) = brint*ux(i,1)
                hgamv(i) = brint*vx(i,1)
           endif
       endif
       prpbl(i) = 1.0
       dthx  = max(thx(i,k+1)-thx(i,k),tmin)
       dqx   = min(qx(i,k+1,1)-qx(i,k,1),0.0)
       hfxpbl(i) = we(i)*dthx
       qfxpbl(i) = we(i)*dqx
!
       dux = ux(i,k+1)-ux(i,k)
       dvx = vx(i,k+1)-vx(i,k)
       if(dux.gt.tmin) then
         ufxpbl(i) = max(prpbl(i)*we(i)*dux,-ust(i)*ust(i))
       elseif(dux.lt.-tmin) then
         ufxpbl(i) = min(prpbl(i)*we(i)*dux,ust(i)*ust(i))
       else
         ufxpbl(i) = 0.0
       endif
       if(dvx.gt.tmin) then
         vfxpbl(i) = max(prpbl(i)*we(i)*dvx,-ust(i)*ust(i))
       elseif(dvx.lt.-tmin) then
         vfxpbl(i) = min(prpbl(i)*we(i)*dvx,ust(i)*ust(i))
       else
         vfxpbl(i) = 0.0
       endif
       delb  = govrth(i)*d3*hpbl(i)
       delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
     endif
   enddo
!
   do k = 1,klpbl
     do i = 1,im
       if(pblflg(i).and.k.ge.kpbl(i))then
         entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2.
       else
         entfac(i,k) = 1.e30
       endif
     enddo
   enddo
!
!     compute diffusion coefficients below pbl
!
   do k = 1,klpbl
     do i = 1,im
       if(k.lt.kpbl(i)) then
         zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
         zfacent(i,k) = (1.-zfac(i,k))**3.
         wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
         wscalek2(i,k) = (phifac*karman*wstar3_2(i)*(zfac(i,k)))**h1
         if(sfcflg(i)) then
           prfac = conpr
           prfac2 = 15.9*(wstar3(i)+wstar3_2(i))/ust3(i)/(1.+4.*karman*(wstar3(i)+wstar3_2(i))/ust3(i))
           prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
         else
           prfac = 0.
           prfac2 = 0.
           prnumfac = 0.
           phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i)
           wscalek(i,k) = ust(i)/phim8z
           wscalek(i,k) = max(wscalek(i,k),0.001)
         endif
         prnum0 = (phih(i)/phim(i)+prfac)
         prnum0 = max(min(prnum0,prmax),prmin)
           xkzm(i,k) = wscalek(i,k) *karman*    zq(i,k+1)      *    zfac(i,k)**pfac+ &
                       wscalek2(i,k)*karman*(hpbl(i)-zq(i,k+1))*(1-zfac(i,k))**pfac
         !Do not include xkzm at kpbl-1 since it changes entrainment
         if (k.eq.kpbl(i)-1.and.cloudflg(i).and.we(i).lt.0.0) then
           xkzm(i,k) = 0.0
         endif
         prnum =  1. + (prnum0-1.)*exp(prnumfac)
         xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
         prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
         prnum =  1. + (prnum0-1.)*exp(prnumfac)
         xkzh(i,k) = xkzm(i,k)/prnum
         xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
         xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
         xkzq(i,k) = xkzq(i,k)+xkzoh(i,k)
         xkzm(i,k) = min(xkzm(i,k),xkzmax)
         xkzh(i,k) = min(xkzh(i,k),xkzmax)
         xkzq(i,k) = min(xkzq(i,k),xkzmax)
       endif
     enddo
   enddo
!
!     compute diffusion coefficients over pbl (free atmosphere)
!
   do k = 1,km-1
     do i = 1,im
       if(k.ge.kpbl(i)) then
         ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k))                         &
              +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k)))                        &
              /(dza(i,k+1)*dza(i,k+1))+1.e-9
         govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
         ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
         if(imvdif.eq.1.and.ntcw.ge.2.and.ntiw.ge.2)then
           if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3.and.(qx(i           &
             ,k+1,ntcw)+qx(i,k+1,ntiw)).gt.0.01e-3)then
!      in cloud
             qmean = 0.5*(qx(i,k,1)+qx(i,k+1,1))
             tmean = 0.5*(tx(i,k)+tx(i,k+1))
             alph  = xlv*qmean/rd/tmean
             chi   = xlv*xlv*qmean/cp/rv/tmean/tmean
             ri    = (1.+alph)*(ri-g*g/ss/tmean/cp*((chi-alph)/(1.+chi)))
           endif
         endif
         zk = karman*zq(i,k+1)
         rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
         rlamdz = min(dza(i,k+1),rlamdz)
         rl2 = (zk*rlamdz/(rlamdz+zk))**2
         dk = rl2*sqrt(ss)
         if(ri.lt.0.)then
! unstable regime
           ri = max(ri, rimin)
           sri = sqrt(-ri)
           xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
           xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
         else
! stable regime
           xkzh(i,k) = dk/(1+5.*ri)**2
           prnum = 1.0+2.1*ri
           prnum = min(prnum,prmax)
           xkzm(i,k) = xkzh(i,k)*prnum
         endif
!
         xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
         xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
         xkzm(i,k) = min(xkzm(i,k),xkzmax)
         xkzh(i,k) = min(xkzh(i,k),xkzmax)
         xkzml(i,k) = xkzm(i,k)
         xkzhl(i,k) = xkzh(i,k)
       endif
     enddo
   enddo
!
!     compute tridiagonal matrix elements for heat
!
   do k = 1,km
     do i = 1,im
       au(i,k) = 0.
       al(i,k) = 0.
       ad(i,k) = 0.
       f1(i,k) = 0.
     enddo
   enddo
!
   do i = 1,im
     ad(i,1) = 1.
     f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2
   enddo
!
   do k = 1,km-1
     do i = 1,im
       dtodsd = dt2/del(i,k)
       dtodsu = dt2/del(i,k+1)
       dsig   = p2d(i,k)-p2d(i,k+1)
       rdz    = 1./dza(i,k+1)
       tem1   = dsig*xkzh(i,k)*rdz
       if(pblflg(i).and.k.lt.kpbl(i)) then
         dsdzt = tem1*(-hgamt(i)/hpbl(i)-hfxpbl(i)*zfacent(i,k)/xkzh(i,k))
         f1(i,k)   = f1(i,k)+dtodsd*dsdzt
         f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
       elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
         xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
         xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
         xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k))
         xkzh(i,k) = min(xkzh(i,k),xkzmax)
         f1(i,k+1) = thx(i,k+1)-300.
       else
         f1(i,k+1) = thx(i,k+1)-300.
       endif
       tem1   = dsig*xkzh(i,k)*rdz
       dsdz2     = tem1*rdz
       au(i,k)   = -dtodsd*dsdz2
       al(i,k)   = -dtodsu*dsdz2
       ad(i,k)   = ad(i,k)-au(i,k)
       ad(i,k+1) = 1.-al(i,k)
     enddo
   enddo
!
! copies here to avoid duplicate input args for tridin
!
   do k = 1,km
     do i = 1,im
       cu(i,k) = au(i,k)
       r1(i,k) = f1(i,k)
     enddo
   enddo
!
   call tridin_ysu(al,ad,cu,r1,au,f1,im,km,1)
!
!     recover tendencies of heat
!
   do k = km,1,-1
     do i = 1,im
       ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
       ttnp(i,k) = ttnp(i,k)+ttend
       dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)
     enddo
   enddo
   if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend) then
     idtend = dtidx(index_of_temperature,index_of_process_pbl)
     if(idtend>=1) then
       dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-thx+300.)*rdt*pi2d
     endif
   endif
!
!     compute tridiagonal matrix elements for moisture, clouds, and gases
!
   do k = 1,km
     do i = 1,im
       au(i,k) = 0.
       al(i,k) = 0.
       ad(i,k) = 0.
     enddo
   enddo
!
   do ic = 1,ndiff
     do i = 1,im
       do k = 1,km
         f3(i,k,ic) = 0.
       enddo
     enddo
   enddo
!
   do i = 1,im
     ad(i,1) = 1.
     f3(i,1,1) = qx(i,1,1)+qfx(i)*g/del(i,1)*dt2
   enddo
!
   if(ndiff.ge.2) then
     do ic = 2,ndiff
       do i = 1,im
         f3(i,1,ic) = qx(i,1,ic)
       enddo
     enddo
   endif
!
   do k = 1,km-1
     do i = 1,im
       if(k.ge.kpbl(i)) then
         xkzq(i,k) = xkzh(i,k)
       endif
     enddo
   enddo
!
   do k = 1,km-1
     do i = 1,im
       dtodsd = dt2/del(i,k)
       dtodsu = dt2/del(i,k+1)
       dsig   = p2d(i,k)-p2d(i,k+1)
       rdz    = 1./dza(i,k+1)
       tem1   = dsig*xkzq(i,k)*rdz
       if(pblflg(i).and.k.lt.kpbl(i)) then
         dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k))
         f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
         f3(i,k+1,1) = qx(i,k+1,1)-dtodsu*dsdzq
       elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
         xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
         xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
         xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k))
         xkzq(i,k) = min(xkzq(i,k),xkzmax)
         f3(i,k+1,1) = qx(i,k+1,1)
       else
         f3(i,k+1,1) = qx(i,k+1,1)
       endif
       tem1   = dsig*xkzq(i,k)*rdz
       dsdz2     = tem1*rdz
       au(i,k)   = -dtodsd*dsdz2
       al(i,k)   = -dtodsu*dsdz2
       ad(i,k)   = ad(i,k)-au(i,k)
       ad(i,k+1) = 1.-al(i,k)
     enddo
   enddo
!
   if(ndiff.ge.2) then
     do ic = 2,ndiff
       do k = 1,km-1
         do i = 1,im
           f3(i,k+1,ic) = qx(i,k+1,ic)
         enddo
       enddo
     enddo
   endif
!
! copies here to avoid duplicate input args for tridin
!
   do k = 1,km
     do i = 1,im
       cu(i,k) = au(i,k)
     enddo
   enddo
!
   do ic = 1,ndiff
     do k = 1,km
       do i = 1,im
         r3(i,k,ic) = f3(i,k,ic)
       enddo
     enddo
   enddo
!
!     solve tridiagonal problem for moisture, clouds, and gases
!
   call tridin_ysu(al,ad,cu,r3,au,f3,im,km,ndiff)
!
!     recover tendencies of heat and moisture
!
   do k = km,1,-1
     do i = 1,im
       qtend = (f3(i,k,1)-qx(i,k,1))*rdt
       qtnp(i,k,1) = qtnp(i,k,1)+qtend
       dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k)
     enddo
   enddo
   if(lssav .and. ldiag3d .and. qdiag3d .and. .not. flag_for_pbl_generic_tend) then
     idtend = dtidx(ntqv+100,index_of_process_pbl)
     if(idtend>=1) then
       dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f3(:,:,1)-qx(:,:,1))*rdt
     endif
   endif
!
   if(ndiff.ge.2) then
     do ic = 2,ndiff
       do k = km,1,-1
         do i = 1,im
           qtend = (f3(i,k,ic)-qx(i,k,ic))*rdt
           qtnp(i,k,ic) = qtnp(i,k,ic)+qtend
         enddo
       enddo
     enddo
     if(lssav .and. ldiag3d .and. ntoz>0 .and. qdiag3d .and.         &
  &               .not. flag_for_pbl_generic_tend) then
       idtend = dtidx(100+ntoz,index_of_process_pbl)
       if(idtend>=1) then
         dtend(:,:,idtend) = dtend(:,:,idtend) + f3(:,:,ntoz)-qx(:,:,ntoz)
       endif
     endif
   endif
!
!     compute tridiagonal matrix elements for momentum
!
   do i = 1,im
     do k = 1,km
       au(i,k) = 0.
       al(i,k) = 0.
       ad(i,k) = 0.
       f1(i,k) = 0.
       f2(i,k) = 0.
     enddo
   enddo
!
   do i = 1,im
     ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2                  &
        *(wspd1(i)/wspd(i))**2
     f1(i,1) = ux(i,1)+uox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i)
     f2(i,1) = vx(i,1)+vox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i)
   enddo
!
   do k = 1,km-1
     do i = 1,im
       dtodsd = dt2/del(i,k)
       dtodsu = dt2/del(i,k+1)
       dsig   = p2d(i,k)-p2d(i,k+1)
       rdz    = 1./dza(i,k+1)
       tem1   = dsig*xkzm(i,k)*rdz
       if(pblflg(i).and.k.lt.kpbl(i))then
         dsdzu     = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
         dsdzv     = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
         f1(i,k)   = f1(i,k)+dtodsd*dsdzu
         f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
         f2(i,k)   = f2(i,k)+dtodsd*dsdzv
         f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
       elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
         xkzm(i,k) = prpbl(i)*xkzh(i,k)
         xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
         xkzm(i,k) = max(xkzm(i,k),xkzom(i,k))
         xkzm(i,k) = min(xkzm(i,k),xkzmax)
         f1(i,k+1) = ux(i,k+1)
         f2(i,k+1) = vx(i,k+1)
       else
         f1(i,k+1) = ux(i,k+1)
         f2(i,k+1) = vx(i,k+1)
       endif
       tem1   = dsig*xkzm(i,k)*rdz
       dsdz2     = tem1*rdz
       au(i,k)   = -dtodsd*dsdz2
       al(i,k)   = -dtodsu*dsdz2
       ad(i,k)   = ad(i,k)-au(i,k)
       ad(i,k+1) = 1.-al(i,k)
     enddo
   enddo
!
! copies here to avoid duplicate input args for tridin
!
   do k = 1,km
     do i = 1,im
       cu(i,k) = au(i,k)
       r1(i,k) = f1(i,k)
       r2(i,k) = f2(i,k)
     enddo
   enddo
!
!     solve tridiagonal problem for momentum
!
   call tridi1n(al,ad,cu,r1,r2,au,f1,f2,im,km,1)
!
!     recover tendencies of momentum
!
   do k = km,1,-1
     do i = 1,im
       utend = (f1(i,k)-ux(i,k))*rdt
       vtend = (f2(i,k)-vx(i,k))*rdt
       utnp(i,k) = utnp(i,k)+utend
       vtnp(i,k) = vtnp(i,k)+vtend
       dusfc(i) = dusfc(i) + utend*conwrc*del(i,k)
       dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k)
     enddo
   enddo
   if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend) then
     idtend = dtidx(index_of_x_wind,index_of_process_pbl)
     if(idtend>=1) then
       dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-ux)*rdt
     endif

     idtend = dtidx(index_of_y_wind,index_of_process_pbl)
     if(idtend>=1) then
       dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f2-vx)*rdt
     endif
   endif
!
!---- end of vertical diffusion
!
   do i = 1,im
     kpbl1d(i) = kpbl(i)
   enddo
!
!
   end subroutine ysuvdif_run
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
   subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,im,km,nt)
   use machine , only : kind_phys
!-------------------------------------------------------------------------------
   implicit none
!-------------------------------------------------------------------------------
!
   integer, intent(in )      ::     im, km, nt
!
   real(kind=kind_phys), dimension( im, 2:km+1 )                                   , &
         intent(in   )  ::                                                 cl
!
   real(kind=kind_phys), dimension( im, km )                                       , &
         intent(in   )  ::                                                 cm, &
                                                                           r1
   real(kind=kind_phys), dimension( im, km,nt )                                    , &
         intent(in   )  ::                                                 r2
!
   real(kind=kind_phys), dimension( im, km )                                       , &
         intent(inout)  ::                                                 au, &
                                                                           cu, &
                                                                           f1
   real(kind=kind_phys), dimension( im, km,nt )                                    , &
         intent(inout)  ::                                                 f2
!
   real(kind=kind_phys)    :: fk
   integer :: i,k,l,n,it
!
!-------------------------------------------------------------------------------
!
   l = im
   n = km
!
   do i = 1,l
     fk = 1./cm(i,1)
     au(i,1) = fk*cu(i,1)
     f1(i,1) = fk*r1(i,1)
   enddo
!
   do it = 1,nt
     do i = 1,l
       fk = 1./cm(i,1)
       f2(i,1,it) = fk*r2(i,1,it)
     enddo
   enddo
!
   do k = 2,n-1
     do i = 1,l
       fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
       au(i,k) = fk*cu(i,k)
       f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1))
     enddo
   enddo
!
   do it = 1,nt
     do k = 2,n-1
       do i = 1,l
         fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
         f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
       enddo
     enddo
   enddo
!
   do i = 1,l
     fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
     f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1))
   enddo
!
   do it = 1,nt
     do i = 1,l
       fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
       f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
     enddo
   enddo
!
   do k = n-1,1,-1
     do i = 1,l
       f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1)
     enddo
   enddo
!
   do it = 1,nt
     do k = n-1,1,-1
       do i = 1,l
         f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
       enddo
     enddo
   enddo
!
   end subroutine tridi1n
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
   subroutine tridin_ysu(cl,cm,cu,r2,au,f2,im,km,nt)
   use machine , only : kind_phys
!-------------------------------------------------------------------------------
   implicit none
!-------------------------------------------------------------------------------
!
   integer, intent(in )      ::     im, km, nt
!
   real(kind=kind_phys), dimension( im, 2:km+1 )                                   , &
         intent(in   )  ::                                                 cl
!
   real(kind=kind_phys), dimension( im, km )                                       , &
         intent(in   )  ::                                                 cm
   real(kind=kind_phys), dimension( im, km,nt )                                    , &
         intent(in   )  ::                                                 r2
!
   real(kind=kind_phys), dimension( im, km )                                       , &
         intent(inout)  ::                                                 au, &
                                                                           cu
   real(kind=kind_phys), dimension( im, km,nt )                                    , &
         intent(inout)  ::                                                 f2
!
   real(kind=kind_phys)    :: fk
   integer :: i,k,l,n,it
!
!-------------------------------------------------------------------------------
!
   l = im
   n = km
!
   do it = 1,nt
     do i = 1,l
       fk = 1./cm(i,1)
       au(i,1) = fk*cu(i,1)
       f2(i,1,it) = fk*r2(i,1,it)
     enddo
   enddo
!
   do it = 1,nt
     do k = 2,n-1
       do i = 1,l
         fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
         au(i,k) = fk*cu(i,k)
         f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
       enddo
     enddo
   enddo
!
   do it = 1,nt
     do i = 1,l
       fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
       f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
     enddo
   enddo
!
   do it = 1,nt
     do k = n-1,1,-1
       do i = 1,l
         f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
       enddo
     enddo
   enddo
!
   end subroutine tridin_ysu
!-------------------------------------------------------------------------------
end module ysuvdif
!-------------------------------------------------------------------------------
